#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include "../common/general.h"
#include "../common/minilog.h"

#define VERBOSE		INFO | RESULT | WARNING | ERROR | CRITICAL
#define LOGGING		INFO | RESULT | PROCESS
#define COMMENT		"Combined HS with DE"
#define LOGFILENAME	"HSDE.hsde"
#define ALGORITHM	"HSDE"

/********************************************************************************
	Harmony Search Parameters
********************/
#define max_iterations		50000
#define hms					50
#define hmcr				0.95
#define par					0.005
#define fw					.00005
#define F					1.5
#define CR					.0005
#define test_function		g02Rosenbrock
#define precision_digits	15
#define sample_step			1e-10
#define eval_stamp			100000
/*
#define dimension				100
#define lower_bound				-1000.0
#define upper_bound				+1000.0
*/
#include "../benchmark/group02fn.h"

/********************************************************************************
	Harmony Search Global Variables
********************/
struct solution{
	double* vector;
	double fitness;
};
struct solution* hm;
struct solution y;
double precision;
int best=-1;
int eval_counter;
char unique_logfilename[256];
#include "../common/common.h"

/********************************************************************************
	Subordinate Functions
********************/
int solComp(const void* a,const void* b) {
	double pa,pb;
	pa=((struct solution*)a)->fitness;
	pb=((struct solution*)b)->fitness;
	if (pa>pb) return 1;
	else if (pa<pb) return -1;
	return 0;
}
void copy(struct solution from,struct solution *to){
	int i;
	to->fitness=from.fitness;
	for(i=0;i<dimension;i++)
		to->vector[i]=from.vector[i];
}

/********************************************************************************
	Harmony Search Main Functions
********************/
void improvise(int index){
	int i;
	for (i=0;i<dimension;i++)
		hm[index].vector[i]=random3_range();
}

void eval(struct solution *sol){
	sol->fitness=func(sol->vector);
	eval_counter++;
}

int init(){
	int i;
	eval_counter=0;
	build_min_vector();
	precision=1.0/pow(10,precision_digits);
	hm=(struct solution*) malloc(sizeof(struct solution)*hms*2);
	y.vector=(double*) malloc(sizeof(double)*dimension);
	for(i=0;i<hms*2;i++)
		hm[i].vector=(double*) malloc(sizeof(double)*dimension);
	for(i=0;i<hms*2;i++){
		improvise(i);
		eval(&hm[i]);
	}
	qsort(hm,hms*2,sizeof(struct solution),solComp);
}
void harmony_consideration(){
	int i,j;
	for(i=0;i<hms;i++){
		for(j=0;j<dimension;j++){
			if(random3()<hmcr)
				hm[hms+i].vector[j]=hm[(int)(random3()*hms)].vector[j];
			else
				hm[hms+i].vector[j]=random3_range();
		}
	}
}
void pitch_adjustment(){
	int i,j;
	for(i=hms;i<hms*2;i++) {
		for(j=0;j<dimension;j++)
			if (random3()<par) hm[i].vector[j]+=(2*fw*random3()-fw);
/*			if (random1()<par) hm[i].vector[j]+=((2*fw*random2()-fw)+random1()*1e-7+random1()*1e-12);*/
		eval(&hm[i]);
	}
}
void differential_evolution(){
	int i,j;
	int r,a,b,c;
	double p;
	double pmin=1e200;
	for(i=0;i<hms;i++){
		r=randint3(dimension);
		a=randint3(hms);
		b=randint3(hms);
		c=randint3(hms);
		if (best!=-1) a=best;
		else while(a==i) a=randint3(hms);
		while(b==i || b==a) b=randint3(hms);
		while(c==i || c==b || c==a) c=randint3(hms);
		for(j=0;j<dimension;j++){
			if(random3()<CR || j==r){
				y.vector[j]=hm[a].vector[j]+\
				F*(hm[b].vector[j]-hm[c].vector[j]);
				if(y.vector[j]<lower_bound || y.vector[j]>upper_bound)
					y.vector[j]=random3_range();
			}
			else{
				y.vector[j]=hm[i].vector[j];
			}
		}
		eval(&y);
		
		if (y.fitness<hm[i].fitness)
			copy(y,&hm[i]);
		if (hm[i].fitness<pmin){
			best=i;
			pmin=hm[i].fitness;
		}
	}
}

int stopping_criterion() {
	return 0;
}

/********************************************************************************
	Display and Debug Functions
********************/
void show_info(FILE* out){
	fprintf(out,"/*\n%s\n*/\n",COMMENT);
	L(INFO)fprintf(out,"Logging Filename                 : %s\n",unique_logfilename);
	fprintf(out,"Algorithm                        : %s\n",ALGORITHM);
	fprintf(out,"Test Function                    : %s\n",func_name);
	fprintf(out,"Number of Variables              : %d\n",dimension);
	fprintf(out,"Variables Lower Bound            : %.15f\n",lower_bound);
	fprintf(out,"Variables Upper Bound            : %.15f\n",upper_bound);
	fprintf(out,"Variables Precision              : %.15f\n",precision);
	fprintf(out,"Maximum Iterations               : %d\n",max_iterations);
	fprintf(out,"Harmony Memory Size              : %d\n",hms);
	fprintf(out,"Harmony Memory Consideration Rate: %.15f\n",hmcr);
	fprintf(out,"Pitch Adjusting Rate             : %.15f\n",par);
	fprintf(out,"Fret Width                       : %.15f\n",fw);
	fprintf(out,"Crossover Rate (CR)              : %.15f\n",CR);
	fprintf(out,"F                                : %.15f\n",F);
	fprintf(out,"-------------------------------------------------------------\n",hms);
}
void show_memory(){
	int i,j;
	for(i=0;i<hms*2;i++){
		printf("[%02d] f=%.1f -- ",i,hm[i].fitness);
		for(j=0;j<dimension;j++)
			printf("x[%d]=%.1f ",j,hm[i].vector[j]);
		printf("\n");
	}
}

int main(){
	int n=0;
	int i;
	double dist;
	double last_best_fitness=1e200;
	double time_stamp;
	char ufname[100];
	clock_t start;
/*	srand(time(NULL));*/
	mt_init();
	get_unique_filename(LOGFILENAME,unique_logfilename);
	minilog_start(unique_logfilename,VERBOSE,LOGGING);
	init();
	V(INFO)show_info(stdout);
	L(INFO)show_info(LOGFILE);
	start=clock();
	while(n<max_iterations){
		harmony_consideration();
		pitch_adjustment();
		differential_evolution();
		qsort(hm,hms*2,sizeof(struct solution),solComp);
		dist=euclid_dist(hm[0].vector,min_vector);
//		VP(PROCESS,"[%06d] f=%.15f\td=%.15f    \r",n,hm[0].fitness,dist);
		VP(PROCESS,"[%07d] f=%.15f\td=%.15f          \r",eval_counter,hm[0].fitness,dist);
		if(fabs(hm[0].fitness-last_best_fitness)>=sample_step || eval_counter%eval_stamp==0) {
//			LP(PROCESS,"%06d\t%.15f\t%.15f\n",n,hm[0].fitness,dist);
			LP(PROCESS,"%07d\t%.15f\t%.15f\n",eval_counter,hm[0].fitness,dist);
			last_best_fitness=hm[0].fitness;
		}
		n++;
	}
	time_stamp=(double) (clock()-start)/CLOCKS_PER_SEC;
//	LP(PROCESS,"%06d\t%.15f\t%.15f\n",--n,hm[0].fitness,dist);
	LP(PROCESS,"%07d\t%.15f\t%.15f\n",eval_counter,hm[0].fitness,dist);
	VP(PROCESS,"\n");
	LOG(RESULT,"-----------------\n");

	LOG(RESULT,"Best Fitness           : %.15f\n",hm[0].fitness);
	LOG(RESULT,"Closest Distance       : %.15f\n",dist);
	LOG(RESULT,"Time                   : %05.2fs\n",time_stamp);
	LOG(RESULT,"-----------------\n");

	for(i=0;i<dimension;i++){
		LOG(RESULT,"x[%d]= %.15f\n",i,hm[0].vector[i]);
	}
	minilog_stop();
	return 0;
}
